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Abstract This paper derives expressions for the growth rates for the random 
2x2 matrices that result from solutions to the random Hill's equation. The 
parameters that appear in Hill's equation include the forcing strength and 
oscillation frequency The development of the solutions to this periodic 
differential equation can be described by a discrete map, where the matrix 
elements are given by the principal solutions for each cycle. Variations in the 
(?fc,Afc) lead to matrix elements that vary from cycle to cycle. This paper 
presents an analysis of the growth rates including cases where all of the 
cycles are highly unstable, where some cycles are near the stability border, 
and where the map would be stable in the absence of fluctuations. For all 
of these regimes, we provide expressions for the growth rates of the matrices 
that describe the solutions. 

Keywords Hill's equation, Random matrices, Lyapunov exponents 



1 Introduction 

This paper considers the growth rates for Hill's equation with parameters 
that vary from cycle to cycle. In this context, Hill's equation takes the form 

^| + [A fe + ? fe Q(i)]y = 0, (1) 
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where the barrier shape function Q(t) is periodic, so that Q(t + At) = Q(t), 
where At is the period. Here we take At = tt, and the function Q is nor- 
malized so that J T Qdt = 1 . The forcing strength parameters q k are a set 
of independent identically distributed (i.i.d.) random variables that take on 
a new value every cycle (where the index k labels the cycle). The parameters 
Afc, which determine the oscillation frequency in the absence of forcing, also 
vary from cycle to cycle (and are i.i.d.). In principal, the cycle interval At 
could also vary; however, this generalized case can be reduced to the problem 
of equation ([1} through an appropriate re-scaling of the other parameters (see 
Theorem 1 of [AB]). 

Hill's equations [HI] with constant values of the parameters have been 
well studied and arise in a wide variety of applications [MW]. The intro- 
duction of parameters that sample a distribution of values is thus a natural 
generalization of this classic problem. Here we refer to the case with constant 
parameters as the "classical regime" of the general case. 

For this class of periodic differential equations, the transformation that 
maps the coefficients of the principal solutions from one cycle to the next 
takes the form 

hk (h 2 k - I)/ g k 
9k hk 



Mi 



(2) 



where the subscript denotes the cycle. The matrix elements are defined by 
hk = D\ (tt) and gk = j/iC 71 ") f° r the kth cycle, where y\ and y% are the principal 
solutions for that cycle. Note that the matrix has only two independent 
elements rather than four: Since the Wronskian of the original differential 
equation ([lj is unity, the determinant of the matrix map must be unity, and 
this constraint eliminates one of the independent elements. In addition, this 
paper specializes to the case where the periodic functions Q(t) are symmetric 
about the midpoint of the period, so that yi(n) = 2/2(f)> which eliminates a 
second independent element [MW] ; this symmetry applies to the applications 
that motivated this work. 

For transformation matrices M. k of the form ©, the eigenvalues A k can 
be used to classify the matrix types [LR]. The characteristic polynomial has 
the form 

X 2 k - 2h k X k + 1 = 0. (3) 

This equation allows for three classes of eigenvalues A&: For \h k \ > 1, the 
eigenvalues are real and have the same sign, and the transformation matrix 
is hyperbolic symplectic; we denote this regime as classically unstable. When 
\hk\ < lj the eigenvalues are complex and the matrix is elliptic; this regime is 
denoted as classically stable. The remaining possibility is for \h k \ = 1, which 
leads to degenerate eigenvalues equal to either +1 or — 1; these matrices are 
parabolic and are stable under multiplication. 

This paper studies the multiplication of infinite strings of random matri- 
ces of the form ([5]), i.e., the product of N such matrices in the limit N — > oo. 
The problem of finding growth rates for infinite products of matrices with 
random elements was formulated over four decades ago [FU, FK] , where ex- 
istence results were given. We recall the key result here for convenience: 
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For a k x k matrix A with real or complex entries, let \\A\\ denote the Frobe- 
nius norm. 

Theorem (FK): Let X 1 , A 2 , X 3 , . . . form a metrically transitive stationary 
stochastic process with values in the set of k x k matrices. Suppose log + | IX 1 1 1 
exists, where log + i = max(logi, 0), then the limit lim7v_ i . 00 HJ^x*- 1 . . . x 1 \ 
exists. 

Determination of the growth rates are thus carried out in the limit of 
large N, and all probabilistic limits given here are meant almost surely. 

A great deal of subsequent work has studied differential equations of the 
form (TTJ) and the growth rates of the corresponding random matrices [CL, PF, 
LGP]. In spite of this progress, there are relatively few examples that provide 
explicit expressions for the growth rates. The goal of this paper is relatively 
modest: It provides (what we believe to be) new analytic expressions for 
the growth rates of random matrices of the form ([2} . These expressions are 
derived for various regimes of parameter space, as described below. 

The outline of this paper is as follows: Section 2 reviews the astrophysical 
background that led us to this topic. Section 3 considers matrix multipli- 
cation for the case where the solutions are unstable in the classical regime. 
Section 4 develops approximations for this regime and provides some numeri- 
cal verification. Section 5 considers matrix multiplication in the regime where 
the solutions are classically stable. In this case, the transformation matrices 
Mk correspond to elliptical rotations and matrix multiplication is stable in 
the absence of fluctuations; random variations in the matrix elements ren- 
der the solutions unstable. The paper concludes (in Section 6) with a brief 
summary of the results. 



2 Astrophysical Background 

The motivation for considering random Hill's equations arose in studies of 
orbit problems in astrophysics [AK]. When an orbit starts in the principal 
plane of a triaxial, extended mass distribution (such as a dark matter halo), 
the motion is unstable to perturbations in the perpendicular direction. The 
development of the instability is described by a random Hill's equation with 
the form given by equation (fTJ). 

To illustrate this type of behavior, consider an extended mass distribution 
with a density profile of the form 

2 2 2 

p=™ with m 2 = + y + (4) 

where po is a density scale. This form arises in many different astrophysical 
contexts, including dark matter halos, galactic bulges, and young embedded 
star clusters. The density field is thus constant on ellipsoids, where, without 
loss of generality, a > b > c > 0. For this density profile, one can find analytic 
forms for both the gravitational potential and the force terms [AK]. From 
these results, one can determine the orbital motion for a test particle moving 
in the potential resulting from the triaxial density distribution of equation 
(|4|). When the orbit begins in any of the three principal planes, the motion is 
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generally unstable to perturbations in the perpendicular direction [AB, AK]. 
For example, for an orbit initially confined to the x-z plane, the amplitude 
of the y coordinate will (usually) grow exponentially with time. In the limit 
of small \y\ <C 1, the equation of motion for the perpendicular coordinate 
simplifies to the form 

T7r + wf.y = where = , , . (5) 

dt 2 yy y y/ c 2 x 2 + a 2 z 2 + bVx 2 + z 2 V ' 

The time evolution of the coordinates (a;, z) is determined by the orbit in the 
original x-z plane. Since the orbital motion is nearly periodic, the z(t)] 
dependence of to 2 represents a nearly periodic forcing term. The forcing 
strengths, and hence the parameters q k appearing in Hill's equation (JTJ, are 
determined by the inner turning points of the orbit (with appropriate weight- 
ing from the axis parameters [a, b, c]). Since the orbits are usually chaotic, 
the distance of closest approach, and hence the strength qk of the forcing, 
varies from cycle to cycle. The outer turning points of the orbit provide a 
minimum value of lu 2 , which defines the unforced oscillation frequency Afc 
appearing in Hill's equation. As a result, the quantity uj 2 can be written in 
the form 



UJ 
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-„ -A fc + Q fc (t), (6) 

where the index k counts the number of orbit crossings. The shapes of the 
functions Qk are nearly the same, so that one can write Qk = qkQ(t), where 
Q(t) is periodic. The chaotic orbit in the original plane leads to different 
values of Xf. and qk for each crossing. The equation of motion (|5|) for the 
y coordinate thus takes the form of Hill's equation (fTJ), where the period, 
forcing strength, and oscillation frequency vary from cycle to cycle. 



3 Matrix Multiplication for the Classically Unstable Regime 

The goal of this work is to find growth rates for solutions of the differential 
equation ([T]). These growth rates are determined by multiplication of the 
random matrices j\4k (from equation [5]) that connect solutions from cycle 
to cycle. These transformation matrices can also be written in the form 



M k = h k B k 



where Bk 



1 Xk4>k 
l/x k 1 



(7) 



where x k — hk/gk and <fik = 1 — l/^t< By virtue of our assumption on the 
variables (qk, Afe), the matrices A4k form a sequence of i.i.d. matrices. In 
this section, we consider the problem of matrix multiplication with matrices 
of the form ([7]). We specialize to the case where the solutions are unstable 
in the classical regime so that \hk\ > 1 and to the case where Xk > 0. We 
also assume that the hk, Xk, and l/xk have finite means. With the matrices 
written in the form ([7)l. the highly unstable regime considered in [AB] can 
be defined as follows: 
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Definition: Given that solutions to Hill's equation (JTJ) are determined by 
transformation matrices of the form ((7]) , the highly unstable regime is defined 
by setting <fik = 1, This specification thus defines a restricted problem. 

We remark that the above regime applies when the matrix elements 
\hk \ ^> 1, which occurs for forcing strength parameters qk 3> 1 [AB2]. 

The growth rates for Hill's equation ((1} arc determined by the growth 
rates for matrix multiplication of the full set of matrices Aik- For a given 
matrix product, denoted here as M.^ N \ the growth rate 7 is determined by 

7= Jim ilog||MW||, (8) 

where the result is independent of the choice of norm 1 1 • 1 1 . We note that the 
growth rate is called the top or largest Lyapunov exponent. 

Equation ((7|) separates the growth rate for this problem into two parts. 
Let the expectation value of a sequence Xk be denoted by 



fc=l 



1 - 

(X k )= hm -Vl t 



Then the first part 7^ of the growth rate is given by 

1 N 

7h = Jim -5>gMHlog|M- (9) 

N—yoo IS * — ' 
k=l 

We limit our discussion to distributions of the hk for which this limit is finite. 
The remaining part of the growth rate is determined by matrix multiplication 
of the Bk- Note that the original differential equation ([1]) is defined on a time 
interval < t < ir, so that the definition of its growth rate includes a factor 
of 7r [MW], whereas the growth rate for matrix multiplication (O generally 
does not [FK]. Ignoring these normalization issues, this paper focuses on the 
calculation of the growth rates for the matrices Aik and Bk- 

The product of TV matrices of type Bk can be written in the form 



N V V 

All X1A12 
'22 



(l/ Xl )S 2l S, 



(10) 



k=l 

where the first equality defines notation and where 

2 JV-1 2 N-1 

i=i 3=1 



rjiV — ± OH — i. 

1 1 



3=1 Tj 3=1 1 j 

Here, the variables rj are products of ratios of the form 

Xp Jl Xf_t2 • • • X ^ n 

r 3 = " 

XlS-i Xj/<y . . . Xy 



(12) 
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The indices are confined to the range 1 < /x«, < N. The additional factors 
aj , bj , Cj , dj are products of the variables <pj , and can be written in the form 

N 

aj = Y[ 4» P k where pk — or 1 . (13) 
fc=i 



Result 1: For the case where \hk\ > 1 for all cycles, and in the limit of large 
N, the eigenvalue of the product matrix is given by the formula 



X = S U 



0{h 



-2N\ 



(14) 



where each of these quantities should be labeled at the Nth iteration. 
Proof: The characteristic equation of the product matrix of equation (|10p 
takes the form 



A 2 - A(27n + S 2 



S11S22 — ^12^21 — . 



(15) 



The final term is the determinant of the product matrix, and this determinant 
is given by the product of the individual matrices, so that 



N N 
S11S22 - £l2E-2,l — - 4>k) = \\ Tj 



k=l 



- K 

— 1 k 



(16) 



k=l 



Given that \hf.\ > 1 Vfc, this term vanishes in the limit N — > 00. As a result, 
the growing eigenvalue of the characteristic equation (|15p simplifies to the 
form A = Six + ^22- d 

Result 2: The four sums that specify the matrix elements of the product 
matrix are not independent. In particular, for the case where \hk\ > 1 and in 
the limit N — > 00, the ratios of the matrix elements approach the form 



2J\2 2J22 

— — = = constant = / . 

-£11 ^21 



(17) 



Proof: As shown above, the determinant of the product matrix vanishes in 
the limit N — > 00, so that in the limit 



^11^22 — ^12^21 



(18) 



The result implied by the first equality of equation (|17[) follows immediately. 

Further, one can show by direct construction that if the relation of equa- 
tion (|17p holds, then the relation is preserved under matrix multiplication. 
Let the product matrix after N cycles have the form 



q(n) 



St f x i^T 
(l/ Xl )S B f2J B 



(19) 



where / is the constant in equation (|17l) . Then the matrix takes the following 
from after the next cycle: 



B (N+l) = 



S T + {x/xi)4>E B xif(S T + {x/xi)4>S B ) 
{l/xi){S B + (x 1 /x)S T ) f(S B + (x 1 /x)St) 



(20) 
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so that the left-right symmetry relation is conserved. □ 

In the above proof we have adopted notation that is used throughout this 
paper: The subscript '1' denotes the values of the parameters (e.g., x-y) for 
the first cycle in the series. Since the results of this problem can be written in 
terms of this starting value, these initial values play a recurring role. The sub- 
script W denotes the values of the parameters (e.g, xn) appropriate for the 
iVth cycle of the series. In iteration formulae, however, we use unsubscripted 
variables (e.g., x) for the next (N + l)st cycle. 

Result 3: In the highly unstable regime, the ratio of Et to Eb has the form: 



St 



x 

Xl 



(21) 



Proof: From our previous results (see equation [19] of [AB]), the product 
matrix after N cycles has the form given by equation flT9"f with / = 1 (in the 
highly unstable regime). After one additional multiplication, we obtain the 
form given by equation (|20l) with / = 1. We thus find 



E T {N+1) 



J7 X W + {x/ Xl )Z B W _ x_ 

Xi 



E B 



(JV+l) 



E B {N) + {xi/x)E T 



(AT) 



(22) 



For each cycle the ratio xj x\ has a different value, so that no limit is reached 
as N — > oo. However, the ratio at any given finite cycle obeys equation (1211) . 
□ 

To derive an expression for the growth rate for matrix multiplication, we 
first define 

S=E U + E 22 ■ (23) 

As shown in the proof of Result 1, the eigenvalue of the product matrix 
approaches S, as defined above, in the limit N —¥ 00. By construction, the 
iteration formula for S takes the form 



S (N+1) = g(N) 



(x/ Xl )<l>E 21 W + ( Xl / x )E 12 W 



AN) 



(24) 



Using the definition of /, Et, and Eb, this expression can be simplified to 
the form 



S {N+x) = g{N) 



(x/x^Eb^ + (x 1 /x)fE T ^ 



(25) 



Result 4: In the highly unstable regime the iteration formula for the eigen- 
value reduces to the form 



(26) 



This result agrees with that of Theorem 2 from [AB] . 
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Proof: In the highly unstable regime <j> = 1, / = 1, and equation (|2"T]) holds 
for the ratio of Et/Eb- The iteration formula of equation (|2"5j) thus reduces 
to 



g(N+l) = ^(iV) 



1 



(x/zi) + (x N /x) 
1 + xn/x\ 



= S (N) 






X\ + X 




X - 




Xi + XN _ 



(27) 



Since the starting value X\ is fixed, the second factor in square brackets 
approaches unity in the limit N — > oo, i.e., 



JV 



lim | 



fe=l 



= 1 



(28) 



The expression of equation (|2T|) thus reduces to that of equation (|2"r?|) . □ 

Motivated by the result of equation (|2"Tj) for the highly unstable regime, 
we write the ratio of matrix elements for the general case in the form 



(N) 



(N) 



XN 
Xl 



-CUN 



so that 



g(N+l) = S (N) 



1 



(x/xi)(f)+ {x N /x)fa N 

f + Ct N (x N /Xi) 



(29) 



(30) 



where the second equality defines J-jy- The parameter incorporates the 
correction due to the matrices not being in the highly unstable regime. Note 
that / approaches a constant value (from Result 2) and x\ is a constant (by 
definition). The iteration factor can be rewritten in the form 



1 + 



X 2 4> + bdlNXN 



x(b + oinXn) 



where 



b = fxi 



(31) 



Theorem 1: The growth rate for matrix multiplication, with products of 
the general form defined through equation (ITU)) , is given by 
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lim 

iV-voo N 



1 N 



fc=i 



i 



x k 4>k + oth-iXk-i 



Xk{^ + OLk-lXk-l) 



where the ak are determined through the iteration formula 

Xk&k + Xk-lCtk-l 



ak 



x k + snan 



(32) 



(33) 



Proof: Note that existence of the required limit holds by the Theorem of FK. 
Equations ([3TJ1 - I3"TT) show that the growth rate is given by 



7 = lim 

JV->oo N 



., A 1 N 

Elog Tk = lim — >^ log 
5 JV->oc N ^ 



k=\ 



k=\ 



x\4>k + ba k -ix k -i 



x k (b + a k ~\Xk-\) 



(34) 
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where this form is exact, provided that the a k are properly specified. This 
issue is addressed below. To complete the proof, we must also show that the 
growth rate is independent of the value of b, so that we can set b = 1 in 
the above formula. The derivative of the growth rate with respect to the 
parameter b takes the form 

N 

&i v 1 sr- 1 d:Fk r?K\ 
lE = ^N2^T k ^' (35) 

k=l K 

which can be evaluated to take the form 

_ j. 1 ^ (a k -ix k -i) 2 - x 2 k 4> k 



-y 

N ^ 



db n^oo N (& + a k ^ix k -i) [x k {b + a k ^xx k ^i) + x\4> k + ba k -ix k -i] ' 

This expression vanishes in the limit. 

To show that the a k are given by equation f|33[) . we start with the result 
of matrix multiplication from equation (I20j) and use the definition of a k from 
equation (j2"9")l ; these two results imply that 

xi s T {k+1) xi z T (k) + , , 

ak+1 x k+ iZ B ^ x k+1 Z B M + { Xl /x k+1 )Z T M ■ {6<) 

We can then eliminate the factors of St and Sb by again using the definition 
of a k from equation (|29p , and thus obtain 



x\ {x k /xi)a k + (x k+ i/xi)<j) k+ i x k a k + x k+1 4> k+ i 

ak+i = 7——, — -, s = : ■ (38) 

Xk+i 1 + {x k /x k+1 )a k x k+1 + x k a k 

After re-labeling the indices, we obtain equation (l33l) . □ 



4 Approximations for the Classically Unstable Regime 



For classically unstable matrices with \h k \ > 1, Theorem 1 provides an exact 
expression for the growth rate. Since the formulae are complicated, this sec- 
tion presents simpler but approximate expressions for the growth rates for 
the case where (j> k are small (Theorem 2) and where the differences 1 — (j> k 
are small (Theorem 3). We also present two heuristic approximations for the 
growth rates for the general problem. 

Theorem 2: In the regime where the variables 4> k are small, 4> k x k <C 1 Vfc, 
the growth rate for the matrix B k tends in the limit of large N to the form: 

7 = log (l + [{l/x k )M k )] 1/2 ) + O {{x k <j> k )) . (39) 

Proof: We first break up the matrix into two parts so that B k = X + A k , 
where I is the identity matrix and where 



. U X k <p k U 7/fc Mn x 



o x k 4> k 




' 


r\k 


l/x k 




Vk 
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Note that the second equality defines r\k = Xk<pk and yk — 1/xfc. We first 
show (by induction) that repeated multiplications of the matrices Ak lead 
to products with simple forms. The products of even numbers N = 2£ of 
matrices Ak produce diagonal matrices of the form 



JV 



pA 





PP 



(41) 



A ( N) = Ai u) = Y[A k 
fc=i 

where the products Pg are defined by 

l t 

Pi = n fe) (y*-^ and p * = n fe-i) • ( 42 ) 



i=l 



i=l 



Similarly, the product of odd numbers N = 2£ + 1 of matrices Ak produce 
off-diagonal matrices of the form 



n 



A (N) = A (2i + 1) = Y[A k = 



k=l 



Qfm 



Qtvi 







(43) 



where the products Qg are defined analogously to the P(. The product of iV 
matrices Bk can then be written in the form 



N 



fc=i 



^2lVl ^22 



(44) 



Without loss of generality, let N = 2£ be even. Then the matrix elements are 
given by 



£11 



N/2 C£ 

EE( p / 

^=0 j=i 



22 



AT/2 C~ 

EE(^ B 

£=0 j=l 



N/2-1 C? e+1 

^12= E E (Of): 

<?=0 j = l 



(45) 



Ar/2-lC™ +1 

^21= E E (tf ), - 

where is the binomial coefficient and where the subscripts on the Pg and 
Qg denote different realizations of the products. 

The eigenvalue An of the product matrix at the Nth iteration is given 
by its characteristic equation, which has the solution 



A 



N 



- j-Eii + £22 + [(^11 - S22) 2 + 4£ 12 £2iriiyi\ 1/2 } 



(46) 



In the limit of large N, we can make the approximation that En w E22 and 
S\2 ~ £21, so that the expression for the eigenvalue takes the form 

N/2C% N/2-lC? e+1 

A N = E 11 + EMviyi} 1/2 = T,E( p ^3+ E E (Qf)j[myi] 1/2 . 

e=o 3=1 e=o 3=1 

(47) 
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In the limit of large N, all the binomial coefficients are large except for the 
first and last one. We can thus rewrite the above equation in the form 

N/2 N/2-1 

An = Y1 c & ((Pf) + <*)+ E ((Of) + ^) tirnf 2 • (48) 

i=o e=o 

If the realizations of the products (Pe)j were independent, the error terms 
eg would vanish in the limit. However, for a given N, the sums contain C^i 
terms, and C^g > N in general, so all of the terms in the sum cannot be 
independent. We then write the products (-P/) and (Qf) in the form 

{P?)+et = (n i ) t (v j ) l {i + n) t , (49) 

and similarly for {Qf) - This form is exact if one uses the proper expressions 
for the eg. Using this result, the expression for the eigenvalue Am becomes 

N/2 N/2-1 



1=0 1=0 



which takes the form 



(50) 



A N = j^C^( m )^( Vj )^{l + e k f'\ (51) 

fc=0 

If we expand this result, we find that 

A N = 1 + Nivj^iy^il + a) 1 ' 2 + C ? (%)(%)(! + e 2 ) + . . . (52) 

Further, by performing an exact treatment of the first order expansion [AB2] 
we find that e\ = 0. This finding allows us to write the product in the form 



A N = 



l + (v 3 ) 1/2 ( yj ) 1/2 + 0( Vj ) . (53) 



N 



The growth rate thus becomes 



7 = log 



l + fe> 1/2 (%) 1/2 l+Ofe)- (54) 



This last expression is valid provided that rjj <C 1 Vj. □ 

Note that to consistent order, we can replace the limiting form of equation 
(|39p with the equivalent, simpler function 

7 ^[(l/* fe )( ??fe )] 1/2 . (55) 

Figure [1] illustrates how well the approximation of Theorem 2 works. 
For the sake of definitcncss, the variables Xk are log-uniformly distributed 
with log 10 Xk € [—2,2]. The 4>k obey the relation <pk = a</>£fc, where £fc is a 
uniformly distributed random variable over the interval [0,1]. As shown by 
the figure, the limiting form of equation (1391) provides an excellent description 
of the calculated growth rate for sufficiently small 
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Fig. 1 Growth rates for small <pk- The variables <j)k are determined through the 
relation 4>k — ct<t> £k, where is uniformly distributed on [0,1]. The solid curve shows 
the growth rate 7 calculated directly from matrix multiplication as a function of 
the amplitude a^. The dashed curve shows the estimate 72 for the growth rate from 
Theorem 2. The dotted curve shows the difference A'y = 72—7. Note that 7 oc ^/a^ 
whereas A-y oc a$. 



Next we consider the case where the correction factors are close to 
unity. In this case the variables (1 - ^) C 1, and we can expand to leading 
order in (1 — cfik). This procedure leads to the following result: 

Theorem 3: Let 70 be the growth rate for the highly unstable regime where 
<f>k = 1. For small perturbations about this limiting case, the growth rate 
takes the form 7 = 70 — £7, where 

<5 7 = lim if ^-f^l + o ((4(1 - 4>k) 2 )) ■ (56) 
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Proof: We again break up the matrix into two parts 
B k = C k - fife Z with Z 



II 1 




(57) 



where here = Xk{l — 4>k) and Ck is the matrix appropriate for the highly 
unstable regime. Note that Z does not depend on the index k. Here we work 
to first order in the small parameter e^. After N cycles, the product matrix 
takes the form 



N 



N 



4 N) =u^= c i N) -j:^ +o(4)-. 



(58) 



k=l 



k=l 



where the partial product matrices V k are given by 




(59) 



We ignore the case where the Z factors appear on the ends - this effect is 
0(1/ N) and vanishes in the limit. The products of the Ck matrices can be 
written in the form 



C 



(N) 



1 X\ 

l/x N Xi/x N 



where 



Z.j'-p 



N 

IK* 

J=2 



XA-l 



. (60) 



where these results follow from previous work [AB]. As a result, the matrices 
d: 



can be evaluated 



(x k + x k+ i)(x k ~i + x k ) 



The product matrix B k N \ given by equation 
be written in the form 



1 X\ 

l/x N xi/x N 



(x k + Xk + l){x k -l + x k ) 

(61) 

to leading order, can now 



C 



N 



N 



(1 - (t>k)xl 



^ (x k + x k+ i){x k -i + Xk) 



(62) 



The first factor is the product of the matrices for the highly unstable regime. 
Since the second factor is a function (not a matrix) its contribution to the 
growth rate is independent of the first factor and represents a correction to 
the growth rate of the form 



5j 



1 N 
lim — > 



(i - Mx 2 k 



N ^-J [x k + Xk+i)(xk-i + x k ) 



+ 0{e 2 k ), 



(63) 



where the equalities hold to leading order. This correction to the growth rate 
has the form given by equation (I56p . □ 
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Fig. 2 Growth rates for <f>k near unity. The variables <f>k are determined through 
the relation <f>k = 1 — A^^kt where is uniformly distributed on [0,1]. The solid 
curve shows the quantity S'y = 70 — 7, where 7 is the growth rate calculated from 
matrix multiplication and 70 is the growth rate for the highly unstable regime (cj>k 
— 1 Vfc). The dashed curve shows the estimate ($7)3 = (70 — 7)3 for the difference 
in growth rate calculated from Theorem 3. The dotted curve shows the error A — 
(^7)3 — (S7. Note that S'y oc A$ whereas the error term A oc {A^) 2 . 



Figure [5] shows the growth rate for small departures from the highly 
unstable regime. The correction factors are taken to have the form (f>k = 
1— A^k, where is a uniformly distributed random variable over the interval 
[0, 1]. The highly unstable regime corresponds to A^ — > 0. The figure shows 
the growth rate calculated from direct matrix multiplication (solid curve) 
and the approximation from Theorem 3 (dashed curve) plotted as a function 
of the amplitude A^. Both curves plot the difference 70 — 7, where 70 is the 
growth rate for the highly unstable regime (where the <pk = !)• 

Since the general case is quite complicated it is useful to have a good 
working approximation for the case where one is not in one of the two regimes 
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cj>k small or near unity. Toward this end, we first show that the values of ak 
have a limited range: 

Result 5: The variables ak are confined to the range </> m i n < < 1, where 
0min is the minimum value of <f>k- 

Proof: We can rewrite the iteration formula (|33[) for ak in the alternate form 

ak = TTK ' (64) 

where we have defined the composite random variable /3k = ak-iXk~i/xk- 
In the present context, < (3k < oo, and we can show that 

doth 

ni >0 (65) 

for all values of (3k- In the limit /3& — > oo, a/. — > 1, whereas in the limit — > 0, 
ttfe - > 0- Hence <p < < 1 for all cycles. But </> > </> m i n , by definition, so 
that </> min < at < 1. □ 

Approximation 1: As a first heuristic approximation, we replace the full 
iteration expression of equation Q33p for ak with the following simplified form 

= j , 66) 

X + x k 

i.e., we use a k — 1 as an approximation for the previous value [keep in mind 
that x is the value at the (k + l)th cycle]. Using equation (1551) to evaluate 
ak in the iteration formula for J~k, we obtain a working approximation for 
the growth rate. Notice that a k appears in the iteration formula for so 
that we must use equation (IBB")) evaluated at fc rather than k + 1. As a result, 
the iteration factor Th involves the random variables Xk from three cycles, 
or, equivalently (since the Xk are i.i.d.) three separate samplings of the vari- 
ables. We change notation so that Xji, Xj2, Xj 3 denote the three independent 
samplings of the random variables x k - Similarly, let 4>j\ 1 4>j2 denote two in- 
dependent samplings of the <pk ■ The iteration formula for this approximation 
can then be written in the form 

J7=1 Zjl<MZj2 + X j3 ) + X j2 {Xj 2 <l>j2 + x j3 ) 
3 %jl [(Xj2 + X j3 ) + X j2 (x 3 2(t>j2 + X j3 )} 

The growth rate for matrix multiplication can then be approximated by 

1 N 
i=i 

where Fj is given by equation (|67p . As a consistency check, for the restricted 
problem where the (p,j n — 1, the iteration factor /Fj reduces to that appropri- 
ate for the highly unstable regime (see equation [26]). 
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Approximation 2: To derive a second approximation for the growth rate, 
we need a better approximation for the a k . If the values of x k and <p k were 
constant, then the ak would approach a constant value given by 

«fe = 2 {(! ~ %k/xk-i) + [(1 - Xk/xk-i) 2 + 4(x fe /j: fc _i)0 fe ] 1/2 | . (69) 

Even though the Xk and <^>fc are not constant, and the ak vary, we can use 
equation (1691) as an approximation to specify the values of ak appearing in 
the exact formula of equation for the growth rate. After using this form 
to specify the ak, and relabeling the indices, the iteration factor takes the 
form 



X 2 kl (j) k i2xk3 + X k 2 {{Xk3 - Xk2) + [(Xk3 ~ X k 2) 2 + ^X k 2Xkz4>k2\ 

H " r~ ■ 

Xkl (Zx k 3 + X k 2 |(Xfe 3 - X k 2) + [(Xks - X k 2) 2 + 4:X k 2Xk34>k2] 1/2 }^ 



Xkl (2xk3 + X k 2 Uxf-3 - X k2 ) + [(Xfc 3 - X k 2) 2 + ^X k 2X k 34>k2\ 12 

(70) 

In the case <f)j n — 1, the iteration factor of equation (1701) reduces to the 
expression for the highly unstable regime (Result 4). 

Figure [3] shows how well these two approximation schemes work. The (j) k 
variables are chosen from the expression cf>k = l — A^k, where is a random 
variable uniformly sampled from the interval < £fc < 1 and where sets 
the amplitude of the departures of the <j) k from unity. The growth rate is 
shown as a function of the amplitude. 

In [AB], we derived a bound on the difference between the growth rate 
for the general case 7 (considered here) and the growth rate in the highly 
unstable regime 70, i.e., 

7o -7 < ^og fa). (71) 

This bound is shown as the dashed curve in Figure [3J The true growth rates 
fall comfortably between this lower bound and the growth rate for the highly 
unstable regime (where the latter provides an upper bound). 

Thus far, this paper has focused on the regime where the transformation 
matrices are classically unstable. Before considering classically stable matrix 
multiplication in the next section, we note the following result that applies 
at the transition between the two regimes: 

Result 6: Consider the matrix transformation that maps the principal so- 
lutions from one cycle to the next. When the matrix elements g k = yi{^) 
vanish, then the remaining matrix elements are h k = yi {it) = ±1. The trans- 
formation matrix M g o for this case is stable under multiplication. 
(The proof is a simple explicit computation.) 



5 Elliptical Rotations and the Classically Stable Regime 

When the principal solutions hk appearing in the discrete map of equation 
([2]) are less than unity, matrix multiplication is stable for the case of constant 
parameters. In the case of interest, however, the parameters in Hill's equation 
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Fig. 3 Validity of approximations of equation (|68[) and equation (|70[) as a function 
of the deviation of 0^ from unity. The upper solid line shows the growth rate for 
matrix multiplication in the highly unstable regime where cj>k = 1. The lower solid 
curve shows the growth rate for the case where 4>k = 1 — ^4<p£fci where is a 
uniformly distributed random variable < < 1. Th e dotted curve shows the 
estimate for growth rate calculated from equation (|68[) using the same sampling 
of the (f>k va riables; similarly, the dot-dashed curve shows the approximation of 
equation (J7DJ . Notice that both of these approximations are almost identical to the 
actual result. The dashed curve shows the lower limit to the growth rate derived 
in [AB]. 



(HJ) and the matrices ([2]) vary from cycle to cycle. This section considers 
the case where the \hk\ < 1, but vary from cycle to cycle, and show that 
instability results. In this regime, the discrete map takes the form of an 
elliptical rotation matrix [LR] as described below. We thus find the growth 
rates for elliptical rotation matrices for the case where the matrix elements 
vary from cycle to cycle. 
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Definition: An elliptical rotation matrix is defined to be 
£(9;L) 



cost/ 
(1/L)sinf 



-Lsinfl 
cos6> 



(72) 



These matrices have the following properties: 

The product of elliptical rotation matrices with the same value of L produces 
another elliptical rotation matrix, also with the same L, 



£(0 1 ;L)£(0 2 ;L)=£([0 1 +0 2 };L) 



(73) 



As a result, the elliptical rotation matrices form a group. 

For fixed L, matrix multiplication is stable. Specifically, the eigenvalues of 

the product of N matrices (with fixed L) have the form 



exp 



exp[±?6»iv] , 



(74) 



where On is the angle corresponding to the group element produced after iV 
matrix multiplications. 

Result 7: When an individual cycle of Hill's equation is stable, specifically 
when \hk\ < 1, the full transformation matrix Aik takes the form of an 
elliptical rotation. 

Proof: Since \hk\ < 1, we can define an angle 9k such that hk — cos Ok- The 
full matrix Aik given by equation (JJJ then takes the form 



M k = 



COS Ok 

9k 



-(sin 2 Ok)/ 9k 
cos Ok 



cos Ok 
(l/ifc)sin k 



-L k sin k 
cos k 



— £k(0k', Lk) , 

(75) 



where we have defined Lk = (sin Ok)/ 9k- As before, we can factor out the 
cos Ok = hk and write the matrix in the form 



M k = cos k 



1 x k (f>k 
l/x k 1 



cos k B k ■ 



where 



Xk 



Lk / tan Ok and 



tan 2 Ok 



(76) 



(77) 



Equation (|77[) thus specifies the transformation between the random variables 
(xk,4>k) appearing in the original transformation matrix and the random 
variables (9k, Lk) in the corresponding elliptical rotation matrix. Note that 
the values of 4>k are strictly negative in this formulation. Otherwise, the 
matrix Bk has the same form as in equation ([7]). □ 

If we let 7b be the growth rate for matrix Bk, then the growth rate 
for the full matrix Aik takes the form 



7M = IB + 



1 



J im T?E l0 S[ C0s6,fe ] 

N^-oo iV — 
k=l 



(78) 
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The exact growth rate for the matrix Bk (see equation [76]) is given by 
Theorem 1. In particular, equations (|3"2"j) and (|3"3"|) remain valid for negative 
values of the <f>k and can be used to calculate the growth rate. 

Result 8: For an elliptical rotation matrix with constant angle 9 and random 
Lk, the growth rate for matrix multiplication vanishes in the two limits h = 
cos 6 — > and h = cos 9 — > 1 . 

Proof: In the limit ft, — > 1 we have sin = 0, and the elliptical rotation matrix 

becomes the identity matrix. As a result, the growth rate vanishes. 

In the other case where h — > 0, s'm9 = 1, and the matrix takes the form 



-Ok 





1/ifc 







In this case, for even numbers of matrix multiplications, say N 
product matrix takes the form 



iV 



k=l 



R B 



(79) 
2n, the 

(80) 



where the matrix elements are given by the products 



n 

fc=i 



L 



2 k 



L<2k-\ 



and 



P L 



n 

fe=i 



-2k-l 



L 



21, 



(81) 



The eigenvalues of the product matrix are given by A = and A = P B . 
For odd N = 2n + 1, the eigenvalue |A| = (P^P^) 1 / 2 . In either case, in the 
limit of large N, the growth rate for matrix multiplication takes the form 



lim 

JV-s-oo N 



1 N 

vE log 



fc=l 



L 



2 k 



^2k-l 



(logi 2 fe) - (logi 2 fe-i) = 0. 



(82) 



The final equality holds because the Lk are independent. □ 

Elliptical rotation matrices are unstable under multiplication when their 
parameters vary from cycle to cycle: 

Theorem 4: Consider an elliptical rotation matrix with variable angle 9^ 
and symmetric fluctuations of the Lk parameter about its mean value Lq. 
The variations are thus written in the form Lk = Lq(1 + rjk), where the odd 
moments (^ n+1 ) = for all integers n. For small fluctuations |?7fc| < 1, the 
growth 7 rate for matrix multiplication takes the form 

1 1 N 

7 = 9 J™, n E lo § cos2 9k + sin2 dk 



fc=l 



1 



1 



+ o(4) 



(83) 



Proof: We first break up the matrix into two parts so that 



£k = 1 cos 9 k + sin 9 k Z k , 



(84) 
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where I is the identity matrix and where 





l/L k 







The product of N matrices £ k becomes 



(85) 



(86) 



fc=i 



=0fe=l \i=l 



where the subscripts on the products denote different realizations. The prod- 
ucts of even numbers I = In of matrices Z k produce diagonal matrices of the 
form 





Z W =Z i2n) = Yl Z2kZ2k _ 1 = { _ iy 



k=l 







(87) 



where the matrix elements P„ and Pjf are given by equation (|81[) . Similarly, 
the product of odd numbers £ = 2n + 1 of matrices Z k produce off-diagonal 
matrices of the form 



z (£) _ Z (2n+1) 



l[Z 2k+1 Z 2k \z 1 = (-1)" 



.fc=l 



-P^ix 




, (88) 



where the P„ are defined previously. Next we write the expectation values of 
these products in the form 



^)=(n 



L 



2j 



\j=l 



J 1 ' 



1 + mj 

+ V2j-1 



1 



1 + r,j 



(89) 



This expression holds because the odd powers of the r]j vanish in the mean, 
and the samples of the different 77's are independent. 

The eigenvalue An of the product matrix at the Nth iteration can be 
written in terms of its matrix elements, i.e., 



^iV = on + (T 2 2 



(90) 



Without loss of generality, let N — 2K be even. The matrix elements <th 
022 = <y are given by 



k ci 



f1K-2r, 



/ k \*=1 / k 



(91) 



m=0 fe=l 



i=l 



where Cf„ is the binomial coefficient and where we have used equation . 
This expression for a contains the even terms of a binomial expansion. We 
can thus write the eigenvalue in the form 



.1 



JV 



Yl [cos6» fc + ism9 k K 1/2 

k=l 



fe=i 



cos 9 k - i sin 6 k K 1/2 



(92) 
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Next we define 

A k = [cos 2 k + sin 2 6 k n] 1/2 and tan a k = K~ 1/2 tan 9 k . (93) 
The eigenvalue takes the form 

A N = 2 (j[ A k) cos (j2 ak ^) ' (94) 
and the corresponding growth rate becomes 
1 1 N 

7 = - Jim — V log [cos 2 9 k + sin 2 6 k 1l} . (95) 

fc=l 

Using the dehnition of 7Z, we obtain the result of Theorem 4. The order of 
the error term follows by comparing equation (|95[) with the leading order 
expansion [AB2]. □ 

In the regime of small r\ k -C 1, the expression for the growth rate reduces 
to the form 

J=\{sm 2 k ){ri) . (96) 

This section shows that instability does not require a finite threshold for 
the amplitude of the fluctuations in L k . Nonzero amplitude leads to instabil- 
ity with growth rate 7 cx (rj k ). Variations in the original parameters (A^,^) 
of Hill's equation lead to fluctuations in the principal solutions (h kl g k ); fluc- 
tuations in the {h k ,g k ) lead to variations in the L k and hence growth. As 
a result, Hill's equation with random forcing terms is generically unstable. 
One notable exception occurs when the h k = or h k — 1 (Result 8). 



6 Conclusion 

This paper provides expressions for the growth rates for the random 2x2 ma- 
trices that result from solutions to the random Hill's equation ([1]). Theorem 
1 gives an exact expression for the growth rate. Theorems 2 and 3 provide 
approximate growth rates for the regimes where the variables <j) k are small, 
and close to unity, respectively. Additional approximations for are given in 
Section 4. When Hill's equation is classically stable, the discrete map that 
governs the solutions has the form of an elliptical rotation matrix (equ. [72]). 
With fixed elements, such matrices are stable under multiplication; variations 
in the L k parameter lead to instability. For small symmetric fluctuations of 
the length parameter L k , the growth rate is given by Theorem 4. 
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